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This paper describes a new Monte Carlo method based on a novel stochastic potential switching 
algorithm. This algorithm enables the equilibrium properties of a system with potential V to be 
computed using a Monte Carlo simulation for a system with a possibly less complex stochastically 
altered potential V. By proper choices of the stochastic switching and transition probabilities, it is 
shown that detailed balance can be strictly maintained with respect to the original potential V . The 
validity of the method is illustrated with a simple one-dimensional example. The method is then 
generalized to multidimensional systems with any additive potential, providing a framework for the 
design of more efficient algorithms to simulate complex systems. A near-critical Lennard- Jones fluid 
with more than 20000 particles is used to illustrate the method. The new algorithm produced a 
much smaller dynamic scaling exponent compared to the Metropolis method and improved sampling 
efficiency by over an order of magnitude. 
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I. INTRODUCTION 

Simulations of complex molecular systems are gen- 
erally carried out using either the molecular dynamics 
(MD) or the Monte Carlo (MC) method. Each method 
has its own merits. The MD method Q, 0, based on 
the integration of the classical equations of motion of the 
particles, is conceptually the simpler of the two. With 
currently available computer power, MD simulations gen- 
erally cannot be carried out for very long time scales for 
very large systems, making the extraction of true equilib- 
rium properties often difficult. The MC method 0,3, on 
the other hand, relies on stochastic dynamics to generate 
members of the desired ensemble. MC has the ability 
to execute non-physical large-scale transitions that are 
impossible in MD and has the potential to reach equi- 
librium much faster. However, devising these large-scale 
transitions that have reasonable acceptance probabilities 
is not always straightforward. 

To search for ways to enable large-scale transitions to 
be carried out with higher probability in MC simulations, 
one must tackle the core problem, which is the complex- 
ity of the interactions among the particles in the system. 
If there were no interactions among these particles, any 
transition, regardless of its scale, would always be ac- 
cepted. When a move is made in MC, the interactions in- 
volving those particles that are being moved will change. 
In general, the larger the scale of the move and the more 
complicated the interactions are, the larger the change 
in the potential becomes. Consequently, large-scale MC 
moves are very unlikely to be accepted. 

A natural questions arises: Is it possible to reduce the 
complexity of the interactions among the particles, for 
instance, by replacing the actual potential V by a less 
complex potential VI One possibility is proposed in this 
paper. With this method, it is indeed possible to re- 
place the original potential V by an arbitrary V . But 
the procedure has to follow a carefully constructed al- 
gorithm to guarantee that detailed balance with respect 



to the original potential is maintained, so that the cor- 
rect statistics are produced. An idea similar to this has 
been exploited in a number of previously proposed Monte 
Carlo methods, such as J-walking 01 [a , M l, simulated 
tempering 0, [t| , parallel tempering |l(il. lllLll2L IT3. cat- 
alytic tempering . multicanonical J-walking [l5| and 
the approximate potential method 0] . But we will show 
that when generalized to multidimensional systems, the 
present method provides flexibilities and potential advan- 
tages that are not available with these previous methods 
and establishes a theoretical framework for the design of 
possibly more efficient algorithms for simulating complex 
systems. 



II. THE SPS IDEA 

Let x be the configuration of a A^-dimensional system 
and V(x) the potential energy divided by the Boltzmann 
constant fog- The statistical weight of each member of 
the canonical ensemble is exp(— V(x)/T), T being the ab- 
solute temperature. A MC algorithm that generates con- 
figurations consistent with their statistical weights can be 
constructed from any set of transition rules, as long as 
the transition probabilities W between every pair x and 
x' satisfy the detailed balance condition: 



-V(x)/T 



W(x 



x') = e- v ^/ T W(x' 



(1) 



The new MC algorithm we propose proceeds as follows: 

1. First, consider changing the system potential V to 
an arbitrary potential V. This "potential switch- 
ing" decision is carried out with a stochastic switch- 
ing probability 



S(x) 



{AV(x)-AV*)/T 



(2) 



with AV(x) = V(x) - V(x) and AV* is a con- 
stant greater than or equal to the maximum value 
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FIG. 1: Illustration of the SPS procedure. (See text for 
details.) 

of AV(x) over all x. By incorporating AV^* into 
Eqn.©, we ensure that S(x) is between and 1. 
Note that if AV* is very large, the resulting S(x) 
will be small, making switching very infrequent. 
For systems with a V that is bounded from above, 
it is always possible to choose a V which yields a 
finite AV*. Systems with an unbounded V will be 
addressed in Sect.lVl 

2. If the switch is made, the configuration of the sys- 
tem is moved from x to x' with transition proba- 
bility W(x — > x 1 ) chosen to satisfy detailed balance 
on the switched potential V, i.e. 

e-v{*)/Tw{x - x') = e-v(x')/Tw{x' -> x). (3) 

Any W, such as Metropolis 0, may be used here 
as long as it satisfies Eqn. J3J). 

3. If the switch is unsuccessful, the configuration of 
the system is moved from x to x' with transition 
probability W(x — > x') chosen to satisfy detailed 
balance on a pseudopotential 

V(x) = V(x) — T ln[l - S(x)} , (4) 

Similar to step 2, any W may be used here as long 
as it satisfies detailed balance on V. 

After the move from x — > x' is made (accepted or re- 
jected), the cycle is over and the algorithm returns to step 
1. This stochastic potential switching (SPS) idea, with 
the relevant potential switching and subsequent MC tran- 
sition probabilities, is illustrated schematically in Fig. ^ 
Obviously, the SPS algorithm can be used alone in a sim- 
ulation (if the moves are ergodic) or mixed with other MC 
moves. 

With the algorithm defined above, it is easy to prove 
that the composite transition probability: 

W(x -» x 1 ) = S(x)W(x -> x') + [1 - S(x)]W(x -» x'), 

(5) 

when substituted into Eqn.Q, indeed satisfies detailed 
balance with respect to the original potential V. There- 
fore, the MC trajectory generated by this SPS idea will 



produce a sequence of configurations {x} that is consis- 
tent with the canonical ensemble for a system with po- 
tential V(x). It is important to emphasize that the choice 
of V is completely arbitrary — the proof works for all V. 
In a real application, one can exploit this arbitrariness 
to select a V that may be either less complex than the 
original V or less costly to compute. For a system where 
the potential is a sum of additive terms V — Vi, often 
the case for many-particle systems, the switching decision 
can be applied to each Vi separately. This generalization 
will be described in Sect.lVl 

In form described above, the SPS idea is concep- 
tually related to J- walking H, IE 0> parallel temper- 
ing [n], El [3 El, and the "approximate potential" 
method |l6j| . In J-walking, the simulation is stochasti- 
cally switched to a configuration sampled from a higher- 
temperature (T 1 ) simulation of the same potential with 
properly chosen transition probabilities. This is essen- 
tially the same as using a potential that is attenuated 
by a factor T/T' as V in SPS. Similarly, in parallel tem- 
pering, the exchange of replicas between two different 
temperatures is equivalent to having one switched to an 
attenuated potential and the other to a higher potential. 
In the approximate potential method, the simulation is 
switched to an approximate potential, and the new con- 
figuration produced on the approximate potential is ac- 
cepted or reject at the end with a "correction" rate that 
is designed to maintain detailed balance with respect to 
the original potential; whereas in SPS, the switching de- 
cision is made before the move, so that the subsequent 
update on V will always be accepted. Even though SPS 
is conceptually akin to these other methods, we will show 
in Sect. that when SPS is generalized to multidimen- 
sional systems, its offers flexibilities and potential ad- 
vantages that are not currently available in these related 
methods. 



III. EXAMPLE: A ONE-DIMENSION MODEL 

We use a simple one-dimensional example to illustrate 
the basic SPS idea. The model we selected was a har- 
monic potential V(x) = ^x 2 , with x confined to within 

the range [—1,1]. We used five different V(x) — \x n , 
with n = 0, 1, 2, 3 and 4 to demonstrate that the en- 
semble averages were indeed invariant with the choice of 
V and the choice of V is hence arbitrary. For each of 
the five V, we sampled x according to the SPS algorithm 
using simple Metropolis moves on both V and V. n = 2 
is a special case, because for n — 2, V — V, so the switch 
is make with unit probability. The SPS algorithm for 
n = 2 is thus equivalent to the normal Metropolis (i.e. 
non-SPS) algorithm. 

The moments (a;" 1 ) at T = 0.2 for m = 1 to 10 are 
shown in Table |U for the five different V. The results are 
clearly invariant with the choice of V to within statistical 
errors, showing that detailed balance is strictly satisfied 
with respect to V for any V. The switching rate Rs is 
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also shown for each V. Rs is in general lower for those V 
(particularly n = 1 and 3) that are very different from the 
original potential V. Notice also that n — corresponds 
to a flat potential. In this case, if the switch is made, the 
potential is completely turned off. 



IV. RELATIONSHIP TO ISING-TYPE 
CLUSTER ALGORITHMS 

A MC method that is capable of executing large-scale 
moves with high probability was proposed by Swendsen 
and Wang [l^ for the Ising model in 1987. This method 
later led to the discovery of a class of methods now col- 
lectively known as "cluster algorithms" [l^. These clus- 
ter algorithms can be shown to be special cases of the 
SPS algorithm. Whereas these cluster algorithms permit 
large-scale MC moves, they are largely restricted to dis- 
crete models and are not generally applicable to molec- 
ular simulations. The SPS method provides a way to 
transfer the ideas behind these cluster algorithms to con- 
tinuous systems. 

In the original Swendsen- Wang algorithm, the interac- 
tion between each pair of Ising spins Oi and Oj is stochas- 
tically deleted with a probability pd = exp[— Jfcaj — 1)], 
where J is the Ising interaction. When an interaction is 
successfully deleted, it has no effect on the subsequent 
MC move. On the other hand, if an interaction is not 
deleted, it is frozen such that the value of this interac- 
tion is constrained to remain constant in the subsequent 
MC move. After all the interactions have been either 
deleted or frozen, the spins break up into clusters of spins 
having frozen interactions. Each cluster can be flipped 
independently of the others. Swendsen and Wang showed 
that this cluster algorithm satisfies detailed balance and 
it produces much faster equilibration compared to the 
conventional Metropolis algorithm ]vf\ - 

The Ising model has potential V — Ylfij) Vij, where 
Vij — —JoiGj and the sum goes over all nearest- neighbor 
pairs. It is easy to show that the Swendsen- Wang algo- 
rithm can be derived from the SPS algorithm by making 
the special choice Vij = 0. As such, the SPS algorithm 
can be considered as a "generalized" cluster algorithm. 
However, we choose not to use this terminology because 
calling SPS a generalized cluster algorithm would im- 
properly imply some geometric origin. Whereas in dis- 
crete systems such as the Ising model there is an obvi- 
ous geometric interpretation for the SPS algorithm, there 
may not be any in more general continuous systems. 



V. GENERALIZATION TO SYSTEMS WITH 
ADDITIVE POTENTIALS 

While the validity of the SPS algorithm is clear for 
the basic case considered in Sect. |HJ the utility of the 
SPS algorithm in this form is rather limited. There are 
several reasons why this formulation of the SPS algorithm 



may not be very practical. (1) The choice of a good V 
is not obvious. (2) Because S(x) in Eqn. [21 is scaled by 
e" AV "/ T , unless V is close to V everywhere, the switching 
frequency will be in general small. (3) For a V that is 
not bounded from above (as in systems with repulsive 
interactions) , there may not be a way to choose a V that 
keeps AV* finite. 

To make the SPS algorithm more useful, we must first 
generalize it to an additive potential V that can be writ- 
ten as a sum of two or more terms. Any V can be decom- 
posed into an arbitrary sum. Some systems, such as those 
with pairwise interactions, have potentials that break up 
naturally into a sum of terms. In other situations, the 
potential may have two or more distinct parts that are 
responsible for different physical phenomena, such as the 
repulsive and attractive part of a Lennard- Jones poten- 
tial |20|. We will see that the usefulness of the SPS algo- 
rithm is related to how V is decomposed. Our formula- 
tion here is inspired by the ideas of Kandel et al. [2l| who 
have provided a generalization of the Swendsen- Wang al- 
gorithm [18( for discrete-state (Ising) models. 

To illustrate the generalization of the SPS algorithm 
to a continuous system with an additive potential, we 
consider a potential with just two terms V(x) = V\{x) + 
V2(x). Extension to more than two terms is straight- 
forward. We can apply the SPS algorithm in Sect. ITTlto 
each of the terms separately, attempting to switch V\ to a 
new V\ and Vi to another Vi with switching probabilities 

Si(x) = eWW-AVV/T and 52(x) = e (AV 2 (,)-AV 2 *)/T 

respectively. Since the potential terms are additive, the 
switching of V\ and V2 are independent of each other and 
can be performed in any order. 

For two terms in the potential, there are four possible 
outcomes of the switching decision. For each one, the 
next part of the simulation will proceed on a different 
potential: (1) if both V\ and V2 are switched, the new 
potential becomes V\ + V2] (2) if both V\ and V2 are 
not switched, the new potential becomes Vi + V%\ (3) if 
V\ is switched but V2 is not, the new potential becomes 
Vl + V2', (4) if V\ is not switched but V2 is, the new 
potential becomes V\ +V2, where Vi and V2 are defined 
as in Eqn. After the switch is made, the system can be 
moved from configuration x — ► x' on the new potential. 
At the end of the move, the original potential can be 
restored to restart the switching process anew. One can 
easily show that with this algorithm, detailed balance 
with respect to the original potential is strictly obeyed 
along any one of the four pathways. This is a direct 
result of the additivity of V\ and V2. The sum over all 
four pathways therefore also obeys detailed balance. 

In the above, we have considered switching both terms 
in V, but this needs not be. In fact, we can apply the 
switching to an arbitrary subset of terms. For example, 
we may consider switching only Vi to V\ and keeping V2 
"alive" . If the switching is successful, the new potential 
becomes V\ +V2; otherwise, it is tq + V^- This scenario is 
equivalent to using a V2 = V% and thus also satisfies de- 
tailed balance. This strategy may be useful, for example, 
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TABLE I: MC results for the model system V(x) = \x 2 at T = 0.2 using 5 different V = \x n . (x m ) are the m-th moments 
measured by SPS-MC. Rs are the observed switching rate for each V . The uncertainty of the last digit is shown in parentheses. 





n = 


n = 1 


n = 2 


n — 3 


n = 4 


(x 1 ) 


-0.0003(8) 


-0.0009(8) 


0.0008(8) 


0.0004(8) 


0.0010(8) 


(x 2 ) 


0.1698(4) 


-0.1698(4) 


0.1702(4) 


0.1701(4) 


0.1701(4) 


(x 3 ) 


-0.0001(4) 


-0.0005(4) 


-0.0001(4) 


0.0001(4) 


0.0003(4) 


(x*) 


0.0716(3) 


0.0723(3) 


0.0720(3) 


0.0720(3) 


0.0721(3) 


(x 5 ) 


-0.0000(3) 


-0.0004(3) 


-0.0001(3) 


0.0000(3) 


0.0001(3) 


(x 6 ) 


0.0416(2) 


0.0424(2) 


0.0418(2) 


0.0419(2) 


0.0420(2) 


(x r ) 


-0.0000(2) 


-0.0003(2) 


-0.0001(2) 


-0.0000(2) 


0.0000(2) 


(x s ) 


0.0283(2) 


0.0291(2) 


0.0285(2) 


0.0286(2) 


0.0287(2) 


(x 9 ) 


-0.0000(2) 


-0.0002(2) 


-0.0001(2) 


-0.0000(2) 


-0.0000(2) 




0.0211(2) 


0.0218(2) 


0.0213(2) 


0.0214(2) 


0.0215(2) 


AV 


0.5 


1. 


0. 


1. 


0.125 


Rs 


0.151 


0.030 


1.000 


0.020 


0.699 



in the case of a V that is not bounded from above. In this 
case, we can decompose the potential into the repulsive 
(unbounded) and attractive (bounded) parts, but apply 
the switching only to the attractive part. 

It should now be clear why decomposing V into many 
additive terms makes the SPS algorithm more practi- 
cal. In the original formulation of the SPS algorithm 
in Sect.[n] the entire V = ^2 e Ve 1S switched to V. The 
probability for the simultaneous switching of all Vi is 
\\i Si. If the potential contains a large number of terms, 
the total switching probability will be small even if each 
individual Si is close to unity. Therefore, switching the 
entire V is almost impossible, but individual terms in V 
can be switched with a much higher probability. 

Coupled with good physical insights, the additivity of 
V can be exploited to devise efficient SPS algorithms that 
may be more optimal than others. Some systems, such 
as those with pairwise interactions, have a natural de- 
composition for V. This may be used to guide the search 
for an optimal breakup. On the other hand, there may 
be a totally unphysical breakup that affords higher effi- 
ciency. The additivity of V offers immense possibilities. 
In the next section, we will illustrate this using a non- 
trivial many-particle example. 

In fact, crude elements of the basic SPS idea have al- 
ready appeared in one of our recent studies on the Monte 
Carlo simulations of imaginary-time path integrals |22| . 
and these ideas have been proven useful for accelerating 
the sampling of stiff paths. The strategy proposed there 
was later implemented in a large-scale path integral sim- 
ulation of superfluid molecular H2 clusters j^. These 
studies motivated us to refine the crude ideas contained 
in those two papers and formulate the more general the- 
oretical framework for the SPS algorithm that has been 
presented in this section. 



VI. EXAMPLE: A LENNARD-JONES FLUID 
NEAR ITS CRITICAL POINT 

The correlation length of a system diverges near the 
critical point. Small local fluctuations of the system 
at large separations become correlated with each other, 
making Monte Carlo simulations extremely sluggish 0|. 
This so-called "critical slowing-down" problem leads to 
extremely long equilibration time for Monte Carlo simu- 
lations that employ only local updates, such as the con- 
ventional Metropolis algorithm. Nonlocal moves can also 
be made using the Metropolis algorithm, but such moves 
are almost always rejected because of the reason given in 
the last section. 

We will demonstrate how the SPS method can be 
used to deal with the critical slowing-down problem in 
a Lennard- Jones fluid. The simulations were carried 
out in a cubic box with periodic boundary condition. 
A Lennard- Jones potential u(r) = 4e [(er/r) 12 — (er/r) 6 ] 
that is truncated but unshifted at r c = 2.5a was used 
for the calculations. Previously, the critical tempera- 
ture T c and density p c for this system were found to be 
k B T c = 1.1853e and p c a 3 = 0.3197 01 . To check scal- 
ing, we vary the box length L from L/a = 10 to 40, using 
up to 20464 particles to maintain a fixed density. Since 
the heat capacity is expected to diverge as \T — T c \ a [25| . 
the slowing-down problem should be manifested in the 
energy measurement. We compared the scaling of the 
autocorrelation time for the energy estimator with the 
box length L of the SPS method against the Metropolis 
algorithm and found a much smaller dynamical scaling 
exponent for the SPS method. 

To implement the SPS method, we break up the to- 
tal potential V — J2i<j u ( r ij)> where is the distance 
between particles i and j, in two stages. First, we can 
obviously break V up into the individual pair interac- 
tions u(rij). Next, for each of the pair interactions, 
we can further decompose it into its positive and neg- 



ative parts, u = u + + U-, such that u+(r) is everywhere 
zero except for r < a where u + (r) = u(r), and w_ is 
its complement. With this, the total potential becomes 
V = J2i<j u +( r ij) + u -( r ij)- (We have also tried to de- 
compose u according to the WCA prescription [2(j but 
found no major difference in the efficiencies of the two 
breakups.) With this decomposition, u + is a purely re- 
pulsive potential, whereas u_ is bounded from above by 
zero. We apply the SPS algorithm to each of the u_(rjj) 
terms to try to switch it to u_ = but keep all the 
u+(rjj) alive. Since it_(r) < u_ for all r, Au*_ can be 
simply set to 0. 

The SPS algorithm, when applied to the Lennard- 
Jones fluid, proceeds as follows. Starting from the cur- 
rent configuration {r*;}, we attempt to switch off each of 
the U—(rij) one by one. After the switching decision has 
been completed for every u_(ry), the particles now in- 
teract with a modified potential in which some pairs of 
particles interact with u- while the rest have zero attrac- 
tion between them. Since u + have been kept alive, every 
pair of particles also interact through the purely repulsive 
At this point, one can employ any Monte Carlo move 
to update the system on this stochastically modified po- 
tential. One simple possibility is to apply a Metropolis 
algorithm with a local update to the SPS modified po- 
tential, just like on the original potential. But as we 
will see, doing this alone will not improve the dynamical 
characteristics of the sampling. 

The autocorrelation time r in MC pass for the to- 
tal energy measurement is shown in Fig. |3 for different 
box sizes L, comparing the conventional Metropolis algo- 
rithm with a local update applied to the original potential 
(square) against the same Metropolis update applied to 
the SPS modified potential (triangles). In both simula- 
tions, one MC pass is defined as having attempted one 
move for each particle in the system. Not surprisingly, 
the dynamic scaling behaviors of the two are identical to 
each other. Since both simulations are based on the same 
local update method and they both satisfy detailed bal- 
ance, the dynamical characteristics of the two sampling 
methods ought to be the same. To improve sampling 
efficiency, nonlocal moves must be used. 

In a near-critical system, the slowing-down problem is 
related to the divergence of the correlation length. Up- 
date methods that employ only local moves will therefore 
have poor dynamical scaling, since they are unable to ef- 
fect large-scale rearrangements of the system. SPS pro- 
vides a basis for designing possibly more efficient alterna- 
tive update schemes. The system potential can be signif- 
icantly simplified using SPS, enabling large-scale moves 
to be performed with much higher acceptance ratio com- 
pared to the original potential. 

To perform large-scale moves in the near-critical 
Lennard-Jones fluid on the SPS modified potential, we 
employed a simple scheme based on an algorithm origi- 
nally proposed by Dress and Krauth [2(| to treat hard 
sphere systems. We changed the method to suit the 
present situation, and our algorithm is illustrated in 
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FIG. 2: Scaling of the autocorrelation time r of the energy es- 
timator with the box length L in the simulation of a Lennard- 
Jones fluid at its critical point. The dynamical exponents r\ for 
the conventional Metropolis method on the original potential 
(squares) and on the SPS modified potential (triangles) are 
identical and equal approximately 2.8. The dynamical expo- 
nent in the SPS algorithm with the cluster reflection update 
(circles) is 1.3. 

Fig. 01 and proceeds as follows. 

1. For every pair of particles interacting with u_, we 
freeze their distance by placing a rigid "bond" be- 
tween them. 

2. For every pair of particles that have a nonzero u+ 
between them, we consider them to be "overlap- 
ping" and also freeze their distance. 

3. Based on the bonds and overlaps, we break the par- 
ticles up into disjoint clusters. Two clusters are dis- 
joint if there is no bond or overlap between them. 
This constitutes the "background" configuration in 
Fig. Ha). 

4. To move the clusters, a point inside the simula- 
tion box is randomly selected to be the pivot (illus- 
trated by the cross in Fig. Eta)), and all particles 
in the box are reflected across the pivot to obtain 
the "foreground" configuration in Fig.^Jb). 

5. When the foreground is overlaid on the back- 
ground, the clusters from the foreground and back- 
ground form additional overlaps (but no additional 
bonds). The foreground and background positions 
of the same particle are also by default considered 
to be in the same cluster. Overlapping clusters can 
then be broken up into disjoint superclusters as in 
Fig. He). 

6. Since there is no overlap between any particles (ei- 
ther foreground or background) from two disjoint 
superclusters, we can choose randomly to accept ei- 
ther the foreground or background positions inside 
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00 (b) 




FIG. 3: Illustration of the cluster reflection algorithm, (a) 
Particles connected by frozen u~ (4-7 and 9-10) are bonded, 
shown in the figure connected by think lines. Particles that 
have a nonzero u+ between them (3-4, 5-6 and 8-10) are shown 
as overlapping. Bonded and overlapping particles break up 
into disjoint clusters. In this example, there are six disjoint 
clusters. This forms the "background" configuration. The 
cross indicates the position of the pivot, which for this il- 
lustration is near the center of the cell, (b) All particles in 
the background are reflected across the pivot to generate the 
"foreground" configuration shown in grey. The new position 
of each particle i in the foreground is labeled i' . (c) Overlay- 
ing (b) on (a) generates superclusters from overlapping fore- 
ground and background clusters. In addition, foreground and 
background positions of the same particle are in the same 
cluster by default. In this example, there are three disjoint 
superclusters: (1', 1, 11', 11), (5', 5, 6', 6), and a third encom- 
passing the rest, (d) For each supercluster, either all the fore- 
ground or all the background positions are accepted into the 
new configuration. In this example, two superclusters take on 
the background positions and one supercluster takes on the 
foreground positions. 



each supercluster with equal probability. An exam- 
ple of the resulting new configuration is illustrated 
in Fig. Eld). 

This completes one pass in our SPS simulation. At this 
point, the simulation can start over with another cluster 
move from step 1, or we can carry out a Metropolis sweep 
before going back to 1. Notice that since the switching 
is done stochastically, a different cluster structure would 
be generated every time even if the switching is applied 
to the same configuration. 

In is easy to show that the cluster reflection algorithm 
above satisfies detailed balance in a trivial way, because 
the move conserves the nonzero part of the total poten- 



tial of the system by fixing all the bonds and overlaps and 
the reflection clearly produces symmetric transition prob- 
abilities. However, by itself the cluster reflection algo- 
rithm is nonergodic, because particles that have no over- 
lap with each other in the configuration before the move 
will have no overlap either after the move. To have an 
ergodic Monte Carlo simulation, this cluster reflection al- 
gorithm must be mixed with another ergodic move, such 
as Metropolis using a local update. In our simulations, 
we performed one Metropolis move for every 10 cluster 
updates, which adds minimal costs to the CPU time. 

Autocorrelation times r in MC pass for the total en- 
ergy measurement is shown in Fig. [21 for the SPS algo- 
rithm using the cluster reflection update (circles) . For the 
SPS algorithm, one MC pass is defined as having made 
one cluster reflection move plus one-tenth of a Metropo- 
lis move (needed for ergodicity). In actual CPU time, a 
SPS MC pass is about 20% faster than a Metropolis MC 
pass. Near the critical point, the autocorrelation time is 
expected to scale with system size as r ~ L v , and the 
dynamical scaling exponent r/ is a measure of the effi- 
ciency of the MC method. Clearly, the SPS method has 
a much smaller dynamical exponent. In terms of absolute 
efficiency, the SPS algorithm is more than ten times bet- 
ter for the largest simulations considered (L/a — 40 with 
20464 particles), and accounting for CPU time difference, 
the SPS algorithm is actually 13 times better. 

A generalized geometric cluster algorithm that is also 
based on the Dress and Krauth cluster move has also been 
proposed recently by Liu and Luijten |27l |. In their ap- 
proach, the energies of the configuration before and after 
the cluster move have to be computed to determine the 
transition probabilities; whereas in our SPS algorithm, 
only the energy before the move has to be computed in 
order to determine the switching probabilities. In cases 
where the potential is complicated and costly to compute, 
the SPS algorithm here will offer CPU time savings com- 
pared to the method of Liu and Luijten, but it may suffer 
from lower switching rates. 



VII. CONCLUSIONS 

In summary, we have presented a new Monte Carlo 
method that is based on a stochastic potential switch- 
ing algorithm. This new algorithm enables the equilib- 
rium properties of a system with potential V to be com- 
puted using a Monte Carlo simulation for a system with a 
possibly less complex stochastically altered potential V. 
Generalization of this method to systems with additive 
potentials provides for an efficient scheme for simulating 
complex systems. The validity of the method is illus- 
trated with a simple one-dimensional example, and its 
practical utility in alleviating the critical slowing-down 
problem is illustrated with a Lennard- Jones fluid near 
its critical point. 
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